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Graphene superlattices (GSLs), formed by subjecting a monolayer graphene sheet to a periodic 
potential, can be used to engineer band structures and, from there, charge transport properties, 
but these are sensitive to the presence of disorder. The localization behavior of massless 2D Dirac 
particles induced by weak disorder is studied for both scalar-potential and vector-potential GSLs, 
computationally as well as analytically by a weak-disorder expansion. In particular, it is investigated 
how the Lyapunov exponent (inverse localization length) depends on the incidence angle to a ID 
GSL. Delocalization resonances are found for both scalar and vector GSLs. The sharp angular 
dependence of the Lyapunov exponent may be exploited to realize disorder- induced filtering, as 
verified by full 2D numerical wave packet simulations. 

PACS numbers: 72.80.Vp, 71. 23. An, 73.20.Fz, 73.20.Jc, 73.21.Cd 



I. INTRODUCTION 

One fundamental aspect of graphene lies in the linear 
dispersion relation of its low-energy charge carriers (elec- 
trons and holes) around the so-called Dirac points. These 
charge carriers behave as relativistic massless chiral Dirac 
fermions and can be described by a two-dimensional 
(2D) Dirac equation. The linear dispersion relation 
is responsible for many discoveries in recent graphene 
research, ' such as half-integer quantum Hall effect, ''' 
Klein's paradox,''' and Zitterbewegung. Other than to 
graphene, Dirac or Dirac-like equations naturally apply 
to cold atoms, ' trapped ions, semiconductors,"' or 
polaritons.'^^ 

Motivated by the importance of Dirac equations in 
such a wide variety of frontier research areas, we study 
in this work disorder-induced localization^''"'"' of massless 
Dirac particles in random potentials. Though our results 
are presented in the context of disordered graphene su- 
perlattices (GSLs, see below) we expect them to be useful 
for many other settings as well. For example, when disor- 
der is introduced to cold-atom simulations of graphene 
or GSLs, our general treatment can be adapted to study 
the impact of randomness on the transport of Dirac mat- 
ter waves. 

GSL refers to graphene under external periodic 
scalar-""-^ or vector potentials. ^ Because GSLs fur- 
ther tailor the band dispersion relation of graphene, 
they may be used to construct graphene-based quan- 
tum devices. Theoretical studies of GSLs and graphene 
under periodic corrugation'' have been highly fruit- 
ful, with remarkable findings such as electron beam 
supercoUimation and the emergence of extra Dirac 
points. On the experimental side, GSLs with scalar 
potential barriers can be created via the electric field 
effect or chemical doping.'""'- Two-dimensional (2D) 
GSLs with a period as small as 5 nm have been created 
through electron-beam induced deposition of carbon. " 



Also triangular GSLs growing on different metal surfaces 
have been observed. " Besides, nano-ripple arrays are 
generated by chemical vapor deposition (CVD)."'' Vector 
potentials are induced by magnetic fields or physical 
strain, ' so vector GSLs can be realized by mounting 
graphene on a substrate with a periodic array of ferro- 
magnetic strips or a periodically structured substrate.''^' ' 

All these laboratory-produced GSLs cannot be per- 
fectly periodic, due to intrinsic randomness and uncon- 
trollable factors during production. Therefore, a more 
realistic GSL should be modeled by a periodic poten- 
tial plus some weak disorder in potential height, poten- 
tial width, or lattice spacing. This randomness causes 
Anderson localization, which turns conductors into in- 
sulators and is especially severe in low dimensions. ^''^'^" 
Consequently, the focus of our work is on the localization 
behavior of a 2D Dirac particle in weakly disordered ID 
GSLs. In a related work, ' ' localization of Dirac parti- 
cles in ID disordered potentials was studied, but only for 
zero incidence angle 9 (i.e., wave vector of charge carri- 
ers normal to the interface between different GSL layers) 
and without analytical results for the localization length. 
Another closely related theoretical study of disordered 
GSLs '" comprised an analytical discussion of the scatter- 
ing transmission only for sufficiently small 9 and random 
barrier heights. Our work extends all previous results, to 
the best of our knowledge, inasmuch as it covers the ana- 
lytical properties of the localization length for all values 
of 9, for different types of disorder, and for both scalar 
and vector GSLs. 

The paper is organized as follows. In Sec. II, we begin 
by modeling disordered scalar and vector GSLs by ID 
rectangular potential barriers or wells. Using a trans- 
fer matrix formalism, we then derive the weak-disorder 
expansion of the localization length, or equivalently the 
associated Lyapunov exponent. In Sec. Ill we present 
analytical and numerical results for the Lyapunov expo- 
nent of scalar GSLs, as modeled by disordered delta or 



2 



rectangular potentials. It is found that at fixed energy, 
the localization length depends very intricately upon the 
incidence angle 9 of 2D Dirac particles in the graphene 
plane. We also predict and confirm the existence of de- 
localization resonances other than for perpendicular in- 
cidence: along these directions the Lyapunov exponent 
vanishes. Our theoretical predictions are fully supported 
by numerical results, as also reported below. Section IV 
is in parallel with Sec. Ill, but treats GSLs with vector 
potentials. In addition, assisted by a numerical study 
of wave-packet dynamics in Sec. V, we propose to use 
the angular dependence of the localization length to re- 
alize a disorder-based filtering mechanism. Section VI 
concludes. 



II. LOCALIZATION LENGTH IN DISORDERED 
GRAPHENE SUPERLATTICES 

A. Disordered graphene superlattices 

Thanks to their linear dispersion relation, low-energy 
charge carriers near the Dirac points in graphene are well 
described by the 2D massless Dirac Hamiltonian: 



5w. 



n+l 



■^n— 1 



Mn-1 



Xn+l 



FIG. 1. Top: Disordered graphene superlattice (GSL) real- 
ized as a scalar potential, Eq. (2), or vector potential, Eq. (6). 
Deviations from the clean GSL (dotted) can occur via fluctua- 
tions in barrier height SV, barrier width Sw, lattice spacing SI, 
and combinations thereof. Bottom: Right- and left-travelling 
wavefunction amplitudes are mapped from one barrier to the 
other by the transfer matrix, Eq. (11). 



Defining Vn = evpAn/c and assuming An{x — Xn) is of 
the same form as Vn{x — Xn), one deals with the same 
parameters as in the scalar case. The different poten- 
tial nature, however, implies very different localization 
properties, as will become clear in Sec. IV. 



H = VF(T-p + V{x). (1) 

In graphene, vp w lO^m/s is the Fermi velocity; cr = 
(fxiCj/) is the vector of Pauli matrices. We consider a 
graphene superlattice (GSL) of parallel potential barri- 
ers, such that the external potential V{x) depends only 
on X. In the following, we consider both scalar and vector 
superlattices. 

A general scalar superlattice potential can be described 

by 

V{x) = Y,ynix-X„). (2) 

We will consider rectangular potential barriers (or wells) 
as depicted in Fig. 1. A perfect GSL has identical po- 
tential barriers (or wells) of height V and width w, i.e. 
Vn{x) = yifO<a;<u' and elsewhere, at lattice posi- 
tions Xn = nl. Due to unavoidable experimental imper- 
fections, or deliberate introduction of randomness, these 
potential parameters fluctuate from site to site: 

Vn = V + 6Vr,, (3) 
Xn = nl+SXn, (4) 
Wn^W + SWn. (5) 

This randomness can induce localization, as will be dis- 
cussed at length in Sec. III. 

A vector-potential superlattice is defined in terms of 
the matrix-valued potential 

V{x) -ay'^'^A„{x - Xn)- (6) 



B. Transfer-matrix formalism 

Because the potential V{x) is separable, the problem 
of describing the transmission across the lattice is effec- 
tively ID, and the transfer matrix formalism is particu- 
larly suited. 

The scattering of a massless Dirac particle through a 
single square barrier (well) is well understood, for scalar 
as well as vector potentials.^'"'' Since the potential is 
piecewise constant, the solution to the Dirac equation is a 
plane wave, both inside and outside the barrier. Outside 
the barrier, solutions of energy E — shvpk with s = ±1 
and k ^ \k\ = (fc^ -I- ky)^/^ are the Dirac bispinors 

^^ix,y) = e^^'^^^+^''^y(^^^l^,ey (7) 

travelling towards right {+) and left (— ), with fc^; > by 
convention. 

= tan-i ^ (8) 

is the incidence angle, or angle of propagation (outside 
the barrier) with respect to the x-axis. 

In the lattice, the wave function between barriers, 
where V{x) = 0, is a superposition of free right- and 
left-moving components created by repeated elastic re- 
flexion and transmission. It is useful to parametrize 
the wave function on the left side of the nth barrier, 
= linie^o+ *(a^n - e)- as 

*n = i^i (,,\«) e-'^y + i^n e^^«^ (9) 
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Since the free solutions (7) between barriers are fixed, 
scattering cannot mix the two components of the 
bispinor, and it suffices to introduce the two ampHtudes 
ip^ , just as for a scalar wave obeying Schrodinger's equa- 
tion on a ID lattice. These amplitudes are mapped from 
n to n + 1 by the transfer matrix: 



(10) 



with 



— e 
t 



tr. 



-iA„ 



-iA„ 



tr. 



(11) 



Reflection and transmission amplitudes r„ and t„ are 
known functions of barrier parameters {Vn,XmWn\ 
and quantum numbers {kx,ky,s\ or equivalently 
{E, 9, s}^'" . A„ = kx{xn+i — Xn) is the free propagation 
phase between superlattice points in the absence of any 
barriers. The transfer matrix is largely determined by the 
symmetries of the scattering problem. ' Unitarity or 
current conservation implies detM„ — 1 — |r„p + |i„p. 
Thus the total reflection and transmission probabilities 
can be expressed as i?„ = |r„p = sin^ (/>„ and r„ — 
|i„P = cos^ (/)„, and we find it useful to parameterize M„ 
as 



e~*^" tan e~*"" sec < 



(12) 



By construction, the net transfer matrix across N barri- 
ers is the product 



N 



(13) 



Before studying this product for the random matrices M„ 
arising from disorder, we first discuss its implications for 
clean GSLs. 



C. Clean graphene superlattices 

In a clean GSL, all transfer matrices Mn — M are iden- 
tical. In other words, a single transfer matrix contains all 
information about the dispersion relation in the lattice, 
which is the essence of Bloch's theorem. 

If parameters are such that |trM| < 2, the energy E 
lies within the conduction band of the GSL. In this case 
the eigenvalues of the transfer matrix are of the form 
A± = e^'^" with /Lt e R, such that trM = 2 cos /x. The 
transfer phase = K^l across one lattice cell determines 
the Bloch vector K^^ of the extended solution in the x 
direction. In terms of the parametrization Eq. (12) the 
dispersion relation in the clean GSL therefore reads'*'' " 



cos KtI 



cos 



sec ©cos a. 



(14) 



The structure of this dispersion is analogous to that of the 
Kronig-Penny model, from which it differs only in the 
functional dependence of the transfer parameters {(j), a\ 
on the potential parameters {V, Z, w} and {£■, 6*, s}. This 
dependence will be made explicit for the two cases of 
scalar and vector potentials in Sees. Ill and IV, respec- 
tively. 

In the case |trM| > 2, the energy E falls into a band 
gap. The wave cannot propagate, and |trM| — 2 cosh k^;? 
defines the exponential decay rate j = KxI across one 
lattice cell. This characteristic localization exponent 
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InlA, 



(15) 



is determined by the larger one of the two eigenvalues of 
the transfer matrix M, which also defines the Lyapunov 
exponent of the product P/y = , whose larger eigen- 
value grows like expl-jN}. 



D. Disordered graphene superlattices 

The transmission across a disordered lattice is de- 
scribed by the product P/v of random matrices shown 
in Eq. (13). Since the pioneering work of Furstenberg,''" 
it is well known that the larger eigenvalue of such a prod- 
uct grows exponentially with probability one. This im- 
plies that a wave incident on the disordered GSL at bar- 
rier number 1 has an exponentially small probability of 
transmission after barrier number N, which is one of the 
hallmarks of disorder-induced localization. Indeed, at 
the first barrier, the wave splits into reflected and trans- 
mitted components, and so on across the lattice. The 
boundary condition is actually simpler after the last bar- 
rier where there is only the transmitted component, 
but no component is incident from the right. Starting 
with the reverse boundary condition (such as tpt ~ 1 
and ipi = 0) at the left, the product in Eq. (13) predicts 
that the solution grows like 



(P. 



N)ll\ 



exp{27V7}, 



(16) 



which suggests the expected exponential localization. 
The Lyapunov exponent, mathematically defined as 



^ = J™ ;^1ii|(Pa,)ii| 



thus determines the localization length Zioc = l/j- Here 
Tjv G [0, 1] is the net transmission probability after N 
barriers. 

The transmission is a random variable, with a very 
wide probability distribution for long enough samples. In 
the localized regime, its most probable (or typical) value 
differs vastly from its mean. The extinction |lnT(iV)|, 
however, has a probability distribution that converges 
towards a normal distribution, such that its most proba- 
ble value is equal to the mean, and the right hand side of 
Eq. (17) indeed converges to the Lyapunov exponent. '"^'''^ 
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While it is an elementary exercise to multiply random 
matrices and extract the Lyapunov exponent numeri- 
cally, there is no simple, general method of calculating 
the Lyapunov exponent exactly for a given model of dis- 
order with arbitrary energy. Different situations require 
different approaches. In the following, we treat two dif- 
ferent cases that are relevant in the GSL context and 
allow for analytical calculations. 



E. Randomly spaced, identical barriers 

First we consider the very simple case where identi- 
cal barriers are distributed with random positions such 
that the free propagation phase between barriers is uni- 
formly distributed in [0, 27r]. Under an ensemble aver- 
age (■) over these random phases, the extinction llnT/vj 
across TV barriers is found to be additive along the sam- 
ple: In Tat = iVlnTi.''''*''^'- Here, Ti is the single-barrier 
transmission at given energy E and propagation angle 
9. Equation (17) then immediately yields the Lyapunov 
exponent 7 = — ^ In Ti . This result holds as long as the 
phases are random enough to satisfy the assumption of 
a uniform distribution, but no matter how small Ti, i.e. 
how strong the scattering. 

At a given energy, a rectangular barrier becomes per- 
fectly transmitting at certain incident angles, and no- 
tably at perpendicular incidence {6 = 0) for all energies — 
this phenomenon is known as Klein tunneling.^ In these 
cases, Ti = 1 implies of course 7 = and absence of lo- 
calization, because all barriers share the same resonance 
condition. 



F. Weak-disorder expansion 

Although the previous elementary model captures the 
essence of disorder-induced exponential localization, it 
cannot describe the more interesting, and arguably more 
relevant, case of barriers with slightly random width, 
height, and/or spacing. In the following, we adapt the 
weak-disorder expansion of Derrida et al.'''* to our case. 
Here, we describe briefly the steps leading to the main 
result; details can be found in Appendices. 

First, we Taylor-expand 



variance: 



M + e,,M'+^M" + 0{el), 



(18) 



where M is the transfer matrix of the corresponding clean 
GSL, the prime (.)' indicates differentiation with respect 
to the perturbed variable V, w or d), and e„ is the weak 
perturbation (e„ = SVn, Swn or 6dn). We assume that 
the random variables at different sites are independent 
and identically distributed, with zero mean and finite 



lim 



e'^ = lim 

N^oc 



1 ^ 

-T 



1 ^ 

n=l 



> 0. 



(19) 



(20) 



Second, we expand the product (13) to order in the 
eigenbasis of A/, where M = diag(A4., A_). The eigen- 
basis of M can be used to find the Lyapunov exponent 
because the exponential growth rate is independent of 
the representation. The matrix element required in (16) 
then reads, neglecting terms of order e'^. 



(Pjv)ii = A 



N 



1 + a; 



N 



(21) 



The second line involves only fluctuations at different 
sites and gives no contribution after the ensemble av- 
erage (see App. Al for details). Inserting the first line 
into Eq. (17), and further using Eqs. (19) and (20), one 
obtains the disorder-induced Lyapunov exponent 
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I 

^loc 



= ^Re X^HJ[\ - (m(^ 



(22) 



In a third step we perform the diagonalization from M 
to M in order to arrive at an explicit expression as func- 
tion of the system variables (see App. A 2 for details). In 
terms of the parametrization (12) one obtains a relatively 
compact result: 



tan 



sin^ fi 



sma 
sin (h 



-\ 2 



■13' 



'2tan2 



(23) 



Here, sin^ is a function of {a,4>} via the clean disper- 
sion relation (14), which is assumed to be satisfied by a 
propagating solution of energy E (otherwise, this pertur- 
bative result of order is merely a small correction to 
the band-gap extinction of Sec. II D). 

Before we discuss the localization exponent (23) in de- 
tail for scalar and vector potentials (Sees. HI and IV), 
we comment on its limit of validity. Eq. (23) diverges at 
the band edges, where sin /i = 0. It is well known that 
localization at these special points occurs with an anoma- 
lous localization length that differs from the perturbative 
result.'' ^''' ' However, exponential localization in the con- 
duction band away from these special points is very well 
described by Eqs. (22) and (23), as will be checked via 
numerical calculations below. 



III. SCALAR POTENTIAL 

Now we specify the transfer-matrix parametrization in 
Eq. (12) for a single scalar potential barrier, the building 
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block of the scalar GSL, in order to analyze the Lyapunov 
exponent given in Eq. (23). The reflection and transmis- 
sion amplitudes r and t are found by piecing together a 
continuous plane-wave solution across the barrier:*^ 

.2 , 



cos if + IS sm 



£ cos 



ql cos ( 



-se 



'""^-e^" tan 6* 



V sm Lp 



(24) 
(25) 



Here e = El/hvp = s\k\l and v = Vl/hvp are energy and 
barrier height expressed in lattice units. Furthermore, 
if — qw is the phase picked up by the plane wave with 
wavevector q = [{v — e)^l~^ — ky]^^^ in the x-direction 
across the potential barrier. 

In the next Sec. Ill A, we first discuss the limiting case 
of 5-like barriers, which admits simple expressions and 
helps to guide the understanding of the general case, 
tackled in Sec. IIIB. 



A. Amplitude-disordered delta scalar potential 

Consider an amplitude-disordered Dirac-Kronig- 
Penney model, made out of regularly spaced (5-peaks of 
random strength, or 5GSL for short. This description 
is appropriate in the low-energy regime, where barriers 
become very narrow and high, k^w 1 and v 3> e. In 
the limit it; — ^ and w — > cx) at fixed vw/l = if, one has 
ql ^ V, and the refiection and transmission coefficients 
in Eqs. (24) and (25) become 



1 sin ip 

— = cos If -\- is- 

r 
t 



cos 6 ' 
-se*^ tan sin 



(26) 
(27) 



These expressions depend on the barrier parameters only 
via the combination ip = vw/l. Therefore, they cover 
randomness in both barrier width and height. We assume 
a resulting phase-shift distribution with mean if = Tp^, 
and fluctuations en = ifn — f with variance = Sip"^. 

Substituting Eqs. (26) and (27) into Eq. (11) and then 
comparing with Eq. (12), one has 



cos If 



cosd , 



e*^ tan (/) = s sin tan 06 e* 



(28) 
(29) 

Taking the real part of the first relation, we find the clean 
dispersion for this (5GSL, 

s 

cos /i — COS kj;l cos f + - — 7 sin k^l sin ip. (30) 

We recall k^l = [e^ — PkyY^"^ = secosO. This relation 
links the energy e to the Bloch-vector components Kj. — 
[ijl and Ky — ky in the bulk superlattice. The dispersion 
is periodic in the potential strength ip, as discussed in 
detail by Barbier et al.*'*', and it suffices to consider < 
ip < 277. In the following, the implications of disorder are 
assessed. 



1. Analytical Lyapunov exponent 

From Eq. (29), it becomes apparent that /3 — k^l — O is 
independent of f. Consequently, /3' = in Eq. (23). Fur- 
thermore, by combining both relations, we can evaluate 
dip{sina/ sm(j)) such that the weak-disorder Lyapunov 
exponent finally reads 



Sif"^ sin^ krl o „ 
7=^^^tan2 0. 
Z sm /i 



(31) 



This expression differs in a characteristic manner from 
the corresponding result for a massive Schrodinger 
particle:'""' instead of being inversely proportional to the 
energy, it is proportional to tan^ 0. This has (at least) 
three important implications: First, instead of diverging 
as the weak-disorder result of Eq. (31) stays valid 
even at low energy e. Second, for perpendicular incidence 
9 ^{), there is no localization, 7 = 0, as required by chi- 
ral symmetry via Klein tunneling^' '. Third, the overall 
angular dependence as tan^ 9 implies that charge carri- 
ers incident with larger angles are localized quite rapidly. 
Therefore, a random (5GSL can act as a directional filter^ 
with preferential transmission perpendicular to the su- 
perlattice barriers. 

Indeed, in the case = 0, i.e., for a purely ran- 
dom potential without a regular superlattice component, 
Eq. (30) reduces to /i = k^l, and the Lyapunov exponent 



7=^tan^0 
' 2 



(32) 



becomes totally independent of energy e. The simple, 
and sharp angular dependence tan^ 9 realizes a disorder 
filter at larger angles, allowing only particles around per- 
pendicular incidence 9 — to transmit ballistically. 

Of course, a richer angular structure arises via the de- 
pendence on kxl = secos9 and the dispersion relation in 
Eq. (30), so that a more detailed discussion is in or- 
der. It is helpful to distinguish two limiting cases. First, 
for e — T' and thus fc^,/ 0, Eq. (30) always has the 
real (hence propagating) solution /i — K^l = f. There- 
fore, unless f = 0,7r, one has sin fc^^Z/ sin/i — > 0, and 
therefore delocalization (7 0) occurs for all angles 
as e — )■ 0. Second, for large enough |e| > tt, one has 
smkxl = for £cos0„ = wk {n ^ 0), while allowing 
cos ^ = ±cosiy9 7^ 0. Equation (31) then indicates that 
sin kxl = gives angular transmission windows 7 = at 



9n ~ arccos(ri7r/£). 



(33) 



Therefore, if the incident angle of a 2D plane wave is 
then the Dirac particle will be delocalized. 



2. Numerical experiment 

We now turn to numerical experiments in order to 
check these predictions. The localization Lyapunov expo- 
nents are extracted numerically by use of Eq. (17), after 
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FIG. 2. Lyapunov exponent 7, Eq. (17), vs. angle 6 at energy 
£ — El/hvp = 27r in a lattice of 5-barriers with ±5% fluc- 
tuations around the average strength (p = 7r/2. In the grey 
shaded areas in panel a, the angle falls into a band-gap sector, 
and 7 is given by Eq. (15). White areas: conduction sectors, 
with 7 due to disorder given by Eq. (31), shown on a mag- 
nified scale in panels b, c. The numerical data confirms the 
delocalization resonance at 9i — arccos(7r/£) — 60°, Eq. (33). 



assume e 



1. The transfer matrix M„ in Eq. (12) 



then depends on the random variable (/?„ via 



M„ = M((^„) = 



cos ipn - 

s sin ipr 



' IS — ^ 
cos P 

tanSe'^ 



s sin ifin tan 9e 



(34) 

The product of transfer matrices obeys the remarkable 
property 



M((^„)A/(V3„_l) = M{ipn + ipn-l)- 



(35) 



Hence, the net transfer matrix across N barriers is Pat = 



Af($jv) where = J2 
transmission probabihty 



N 



ifn- As a consequence, the 



cos^ 6 



COS^ $ 



N 



sin^ <I>f, 
cos^ 9 



1 - sin 6 cos2 $ 



(36) 



N 



is bounded from below by cos^ 6. So for the resonance 
angles 0„ of Eq. (33), Tjv cannot be an exponentially 
decaying function of N, thus proving 7 = 0. We em- 
phasize that this delocalization is no longer based on a 
weak-disorder expansion. Rather, it is an exact result for 
arbitrary disorder strength. 



first multiplying random matrices according to Eq. (13). 
Unless specified otherwise, we always take N = 1000 ran- 
dom potential barriers and then ensemble-average over 30 
samples to reach negligible statistical error. 

Figure 2 compares the analytical result of Eq. (31) 
with the numerical data, for a varying incident angle 9 
at fixed energy e = 27r. The average lattice strength is 
set to ip = 7r/2, and we allow 5% equiprobable fluctua- 
tions (Jw^/w^ — 0.01/12). In the overview panel a, grey 
shading shows the intervals where e falls into a band gap. 
The incident wave then turns into an evanescent wave, 
whose attenuation is described by 7 of Eq. (15), with 
negligible corrections due to disorder. For the given pa- 
rameters, the band edges are located at angles 6b solving 
cos6'f, = ±sin(£cos6'6), i.e. 9b G {32.7°, 52.9°, 64.6°}. At 
these points, Fig. 2 shows hardly visible spikes, where 
the perturbative result of Eq. (31) is expected to fail.''"''''''' 
Inside the conduction intervals, shown in the magnified 
view of panels b and c, the Lyapunov exponent 7 given by 
Eq. (31) is in excellent agreement with numerical results. 
The delocalization resonance at 9i — 60° is confirmed, 
with 7 vanishing there. 

3. Exact delocalization resonance 

Interestingly, the numerical evidence suggests that the 
delocalization resonance not only holds perturbatively to 
order e^, as predicted by Eq. (31), but instead is an exact 
resonance. So we seek non-perturbative insights by re- 
turning to the transfer matrices. Under the resonance 
condition, kxl = e cos 0„ — nn, the free propagation 
phase is 6**^==' = ±1. Without losing generality, let us 



B. Disordered square scalar potential 

We return to the general case of a rectangular potential 
superlattice, and proceed as previously. With Eqs. (24) 
and (25) used in (11), the comparison with Eq. (12) yields 

/ SI) tyP" \ 

e'" sec = e** I cos — i sin (p — ) , (37) 

V 'f?' / 

e^'^tan</. = e^(^-^)tan0'^. (38) 

sql 

We denote k ~ Ik^ = se cos 9 and ql — [v^ — 2ev + n"^]^^^, 
as well as ip = qw. We have introduced 5 — k^d as the 
phase picked up over the distance d — I — w between 
barriers on average. In terms of these parameters, the 
dispersion relation of the clean GSL reads"'' 

SV 

COS /Lt = cos (5 COS (/J -I , — sin (5 sin (39) 

nql 

1. Lyapunov exponent 

In the disordered case, Eq. (38) fixes (3 = k^d — 6, 
which now depends on the distance d — I — w between 
consecutive barriers, such that 13' in Eq. (23) is finite for 
barriers of variable distance d. By combining Eqs. (37) 
and (38), one finds 

f sin^ tp ^,,1 sin^ 09 9^ 

2 iq^PsiTL^ ) q'^l'^ 



7 




reduced to 



FIG. 3. Lyapunov exponent 7, Eq. (17), vs. angle 6 at energy 
e — El/hvp = 2-K in a lattice of potential barriers with fixed 
width w = 0.51, and 5% fluctuations around average height 
V = Vl/hvp = TT, such that wv/l = n/2 matches the S barrier 
strength Lp used in Fig. 2. In the grey shaded areas in panel 
a, the angle points into a band-gap direction, and 7 is given 
by Eq. (15). White areas: conduction sectors, with 7 due 
to disorder given by Eq. (40), shown on a magnified scale 
in panels b, c. The numerical data for stronger disorder in 
panel c shows that the perturbative delocalization resonance 
at 6^ ~ 51.5° exists only for weak disorder. 



where S" denotes the derivative of 



S = 



ql sin 6 cos ip + k^l cos S sin ip 
V sin (y5 



cos (5 (41) 



with respect to the fluctuating barrier parameter. 

The relation between the Lyapunov exponent 7 and 
the incidence angle 9 established by Eqs. (40) and (41) is 
rather complicated. But the overall factor tan^ 6 guaran- 
tees absence of localization due to Klein tunneling in the 
forward direction, as well as the efficient filtering of large- 
angle components. To verify our analytical prediction, 
we plot in Fig. 3 the Lyapunov exponent 7 as function 
of 9 at fixed e for barriers with randomly varying height, 
together with the data from the numerical solution. The 
agreement is excellent. 



2. Approximate delocalization resonance 

The numerical results show that there exists a delocal- 
ization resonance 7 = also in this case, but at a slightly 
different angle, 9 « 51.5°, compared to the (5-barriers of 
Sec. Ill A. Let us see how this result comes about. With 
/3' = (because in Fig. 3 only the barrier height fluctu- 
ates) in Eq. (40), there are two non-trivial factors that 
can vanish, simp and S". 

First, there is the obvious candidate sin ip — 0, which is 
the single-barrier resonance condition qw ~ rnr. But this 
zero is exactly cancelled by the most singular contribu- 
tion to S' in Eq. (41), namely S" — (3'(/3'sin(5/(wsin^ ip) + 
O [(sintp)"^] . Since for sm(p — the dispersion relation 
Eq. (39) reads cos /i = ± cos S, Eq. (40) can be further 



Sv^ w'^v'^{v — e)^ 



tan^ 9. 



(42) 



This expression could be thought to vanish for v — e, i.e., 
when the energy equals the mean potential height. How- 
ever, V = s implies q = iky, which is impossible because 
it contradicts the initially assumed resonance condition 
qw = mr. 

Therefore, S" = must be responsible for the observed 
delocalization resonance 7 = 6. In general, the equa- 
tion S" = is too complicated to admit an analytical 
solution, but the resonance angles 0„ can be found nu- 
merically. For the present parameters it is the resonance 
angle 9i that is observed in Fig. 3. In contrast to the 
case of (5-barriers, though, this resonance is not exact. 
In Fig. 3c, numerical results for stronger disorder show 
a deviation from 7 = 0, thus indicating the absence of a 
true delocalization resonance. 

We note that for a purely random potential, v — and 
En = Vn, the Lyapunov exponent Eq. (40) reduces to 



~2 kW 



tan 



(43) 



Now the single-barrier resonance condition k^w = wk 
does lead to 7 = 0. Consistently, the limit d — )■ of 
Eq. (42) vanishes. Here, to lowest order in e = 5v, the 
wavevector inside the barrier is kx, and the resonance 
condition can be satisfied everywhere. But it needs to 
be emphasized that also this result holds only for weak 
disorder, and hence the Lyapunov exponent is not abso- 
lutely zero due to higher-order terms of e. 



IV. VECTOR POTENTIAL 

This section parallels the previous one, with results 
pertaining to disordered vector-potential GSLs, as intro- 
duced in Eq. (6). The single-barrier reflection and trans- 
mission amplitudes r and t are ' ' 



1 iwk f ■ ■ usiTi9 + ecos.- 6 , 

- = e™"" cos <y5 - I sm 03 ; , (44) 

t V g cos 6* ; ' V ^ 



t 



e sec( 



u sm ip 



(45) 



Here e = El/hvp = s\k\l and u = eAl/hc are energy 
and barrier height expressed in lattice units. Besides, 
ip — qw is the phase picked up by the plane wave with 
wavevector g = l~^[e^ — {lky — u)'^]^^^ across the potential 
barrier. The variable q differs from the wavevector q in 
the previous scalar potential case. In particular, q can 
be imaginary if u is large, leading to bound states inside 
a barrier. ' 
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A. Amplitude-disordered delta vector potential 

Very narrow and high potentials barriers i.e., kxW <C 1 
and u e, reahze a vector (5GSL. In the hmit w — > 
and M — >■ oo at fixed uw/l — i.p, one has q — >■ iu/l, and 
the reflection and transmission coefficients in Eqs. (44) 
and (45) reduce to 



— cosh Lp ~ i sinh Lp tan 9, 
= e'* sec 9 sinh ip. 



(46) 
(47) 



The fluctuating phase (p — uw jl describes randomness in 
both width w and height u. We assume a distribution 
with mean ^ — , and small fluctuations — ip^^ — ^ 
with variance = bip"^ . 

Substituting Eqs. (46) and (47) into Eq. (11) and com- 
paring with Eq. (12), one has 



e*" sec = (cosh ip + i sinh ip tan 6l)e'''"' , (48) 
- sinh sec 6'e-'^e*'=-'. (49) 



e'^ tan 



The clean dispersion (14) for the vector 5GSL is found 
by taking the real part of the first relation: 

cos /i = cosh Lp cos K — tan 9 sinh ^ sin k. (50) 

where k — k^l = [e^ — Pky]^/^. 

1. Lyapunov exponent 

For regularly spaced potentials, /3' = in Eq. (23). 
Using Eqs. (48) and (49) to evaluate 9^ (sin a/ sint/)), we 
find the weak-disorder Lyapunov exponent 



S(p'^ sin kxl 2 

7 = — sec I 

z sm ^ 



(51) 



This expression ressembles much the scalar (5GSL result 
(31), except that the tan^ factor is replaced by sec^ 6*. 
Therefore, vector (SGSL and scalar (5GSL share (at least) 
one interesting feature: Instead of diverging as e^^ in the 
Schrodinger case, the weak-disorder prediction of the lo- 
calization length stays valid even at low energy e. On the 
other hand, because sec^ = 1 at 6* = 0, there is no rea- 
son to expect delocalized solutions close to perpendicular 
incidence on general grounds. 



2. Absence of delocalization resonances 

For ip = 0, representing a random vector potential with 
zero mean, Eq. (51) reduces to the energy- independent 
expression 7 = ^Sip'^ sec^ 9. The angular dependence 
sec^ 9 differs from Eq. (32) for scalar (5GSL in that local- 
ization stays finite even at perpendicular incidence 9 = 0, 



2.5" 

2 
1.5 

1 < 
0.5 



a 





FIG. 4. Lyapunov exponent 7, Eq. (17), vs. incident angle 6, 
for a vector (5GSL, with average lattice strength ip = 1, energy 
e = 2tt, and disorder modeled by ±5% fluctuations around ip. 
In the grey shaded areas in panel a, the incident angle falls 
into a band-gap direction. The white area is the conduction 
sector, in which the anaytical Lyapunov exponent is given by 
Eq. (51), as shown in Panel b on a larger scale. 



but becomes just as strong at grazing incidence 9 —> 7r/2 
where sec 6* « tan6'. 

For the general situation with (/? ^ 0, we analyze two 
representative cases. Consider first the low-energy limit 
e — >■ and thus k — 0. Then the dispersion Eq. (50) 
reads cos/x = cosh(/3, which requires an imaginary /i 
and hence describes a non-propagating solution inside 
the band gap. As such, the vector (5GSL acts as an insu- 
lator for small e and arbitrary incidence angle 9. 

Next we turn to cases with sufficiently large |e| > tt. If 
£ cos 9n = {n 0), one has sin k = and hence 7 = 0. 
Note that under this condition, |cos/i| = cosh(/?, which 
again implies a non-propagating solution. So here 7 = 
merely indicates that the disorder-induced correction 
to the decay exponent of the evanescent wave is zero. 
Putting all the above considerations together, it appears 
that the localization behavior in vector (5GSLs is not as 
rich as in scalar (5GSLs. 

Figure 4 compares the analytical prediction Eq. (51) 
with numerical results, for varying incidence angle 9 at 
fixed energy e = 2tt. The average lattice strength is 
ip = 1, and disorder is modeled by 5% equiprobable fluc- 
tuations around ip. In the overview panel a, band-gap 
regimes are grey shaded. The band edges lie at the an- 
gles 9b & {42°, 57°}. Exactly at these points, abnormal 
spikes are seen in panel a, signaling the expected failure 
of Eq. (51).'''"'' ' Inside the conduction band, shown on a 
magnified scale in panel b, the agreement between theory 
and numerics is excellent. In particular, no delocalization 
resonance is seen, as analyzed above. 
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B. Disordered square vector potential 

The general case of a disordered rectangular vector po- 
tential is somewhat more complicated and considerably 
richer in physics. Following the same procedure as for 
scalar GSLs, the first step is to connect the transfer- 
matrix parameters a, (3 and (j) to the GSL parameters 
w, u, and I, as well as the Dirac-particle quantum num- 
bers £, 9, s. For that purpose we use Eqs. (44), (45), (If), 
and (f2) to obtain 

,A / + ue sin 9\ . . 
e sec = e cos (/3 -I- i sm (/3 — , (52) 



IqK 



e"^ tan, 



-e^(^-«)^sec( 
Iq 



(53) 



flere k = Ik^ — ecos9, Iq — [e^ — [Iky — u)^]-'^/^, and 
Lp — qw. 5 — k,j.[l — w) is the phase picked up between 
neighboring barriers. In terms of these parameters, the 
dispersion relation of a clean GSL becomes 



cos IX — cos 5 cos 



ue sm f 



Iqn 



sin 5 sin (p. (54) 




1. Lyapunov exponent 



FIG. 5. Lyapunov exponent 7, Eq. (17), vs. incident angle 9, 
at energy e — 2n for a disordered lattice of rectangular vector 
potentials. Disorder is modeled by ±5% fluctuations around 
an average barrier height u = 2, while periodicity I and barrier 
width w — 0.51 are fixed, such that wu/l — 1 matches the 5 
barrier strength 95 used in Fig. 4. The grey shaded area in 
panel a indicates a band gap and the white area indicates 
conducting solutions. Panels b and c show details of the 
conduction sector. An approximate delocalization resonance 
appears around 9 — f8.5°, for both analytical and numerical 
results. The numerical data for stronger disorder in Panel c 
proves the departure from the weak-disorder resonance. 



Using Eqs. (52) and (53), we can apply our general 
result Eq. (23) once again, leading to 

^-IItS^I^T+Z^'^^-^^, (55) 
2 [Pq^sm 1.1 J l''q'' 

where S' denotes the derivative of 

~ sin 5 cos <Z3 + K cos (5 sin 03 e c ■ ^ 

b = : 1 — cosd sm0 (56) 

u sm if K 

with respect to the fluctuating barrier parameter. In 
contrast to the scalar potential with overall tan^ 9 de- 
pendence, the factor sec^ 6 in Eq. (55) does not lead to 
a simple delocalization resonance at perpendicular inci- 
dence, just as for the vector (5GSL of Sec. IV A. 

Our numerical data confirm these predictions, as seen 
in Fig. 5. The statistical fluctuations in Fig. 5b appear 
larger than before because for the present parameters, 
the Lyapunov exponent 7 is extremely small. 

2. Approximate delocalization resonance 

Figure 5 also reveals a delocalization resonance 7 = 
at 6* w f 8.5° , all the more remarkable because no such res- 
onance occurs in the (5-barrier limit of Sec. IV A. In order 
to explain this analytically, we return to Eq. (55). First 
of all, for the amplitude randomness studied in Fig. 5, 
l3' = 0. Then, 7 = at u 7^ imphes sirnp^^O or S' = 0. 

Let us begin by analyzing the case simp = 0, which 
is equivalent to the barrier resonance condition qw = 



{w/l)[e^ — (esin6' — u)"^]^/^ = mr. To leading order in 
f / sin we find |>S"| = \q(p' sin (5/(m sin^ ip) \ from Eq. (56). 
So the sin* ip factors cancel in Eq. (55), which reduces to 

5V? w'^U^{u-lky)'^ 2/1 /r-TN 

This expression vanishes (remember u 7^ 0) for u = Iky, 
which is equivalent to m = e sin 9. Together with the bar- 
rier resonance condition, this fixes e„ = rnrl/w. There- 
fore, resonances should occur whenever 

9n = arcsin(M/e„) (58) 

For the parameters of Fig. 5 {u = 2, e = 27r, and 
w = 1/2), Eq. (58) predicts a resonance at 9i f8.6°, 
in perfect agreement with the data in Fig. 5c. As shown 
by the data for stronger disorder, the delocalization res- 
onance only holds to lowest order of the weak-disorder 
expansion. 

Are there other delocalization resonances caused by 
S" = 0? A direct answer is difficult on account of the 
rather complex expression for S". Numerically, we have 
scanned the values of S' and find that when S' is zero, 
the associated solution falls inside a band gap. This being 
the case, the S' — condition does not produce new delo- 
calization resonances, in marked difference to the scalar 
GSLs studied in Sec. IfIB2. 
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FIG. 6. Contour plot of the probability density, from 2D 
wave-packet dynamics simulations, for various amplitudes 
V = Vl/hvF of clean (upper row) and amplitude-disordered 
(lower row) scalar GSL potentials V{x), Eq. (2). The time 
evolution samples all incidence angles 6 at once, starting 
with an isotropic wave packet centered around energy e = 
El/hvF = 2n (see text). Comparison between panels a-b, 
and c-d demonstrates disorder-induced filtering: since wave- 
packet components at larger angle 9 have a shorter localiza- 
tion length, they cannot contribute to propagation in x di- 
rection, and the transmitted part of the wave packet appears 
more focused. The GSL potential in panels e-f is sufficiently 
strong to induce the wave-packet collimation that accompa- 
nies the emergence of new Dirac cones. Panel f shows that 
disorder has rather little effect on collimation. 



V. WAVE PACKET DYNAMICS: 
DISORDER-INDUCED FILTERING 

Our analytical results have revealed an interesting 
functional dependence of the localization length upon the 
incident angle of charge carriers. In particular, the Lya- 
punov exponent 7 = l/koc of a scalar GSL is proportional 
to tan^ 9. This factor indicates a strong angular depen- 
dence of disordered-induced localization: the localization 
length diverges for small 9 and quickly decreases as 9 in- 
creases. Certainly, for 9 too close to 9 = tt/2, an infinite 
Lyapunov exponent or vanishing localization length is 
an artifact of weak-disorder perturbation theory. With 
this clarified, it is nevertheless clear that scattering waves 
with larger 9 tend to be much more localized than those 
with small 9. And wave components with localization 
length shorter than the GSL sample will not contribute 
to the conductance. This realizes a filtering effect due to 
disorder. The main goal of the present, comparatively 
short section is to confirm this effect by a direct dynami- 
cal simulation of wave-packet transmission across a scalar 
GSL, both with and without disorder. 



Figure 6 shows the result of a numerical solution of the 
time-dependent Dirac equation with Hamiltonian (1) and 
a scalar GSL potential, Eq. (2), with symmetric barrier 
width w — 0.51 filling the half-space x > 101. In order 
to sample all incidence angles at once, we choose as ini- 
tial condition an isotropic wave packet with momentum 
components ^'(p) oc exp{ — (|p| — po)^/(2Ap^)} centered 
on the radial value po = 2Tih/l with spread Ap = 0.2h/l] 
the wave packet's central energy therefore is £ = 27r in 
lattice units. In Fig. 6, we plot the probability density 
at time t = 701 /vp\ in some cases, a substantial part of 
the wave packet is reflected into the half-space x < 101 
(not shown). The upper row shows the results for clean 
GSLs of different strengths, whereas the lower row shows 
the results for a single realization of the corresponding 
disordered GSLs with fluctuating potential heights. 

Panels a and b compare a pristine graphene sheet 
to a purely amplitude-disordered scalar GSL with zero 
mean potential strength and equiprobable fluctuations 
6v e [—1,1]. Whereas the clean substrate allows for 
isotropic propagation, in the disordered GSL the larger- 
angle components are localized more strongly, as ex- 
pressed by the tan^ 0-behavior of the Lyapunov expo- 
nent, Eq. (43). Consequently, the propagating part of 
the wave packet is concentrated around the forward di- 
rection 9 — 0^ thus supporting our filtering conjecture 
above. 

Panels c and d compare again the clean and disordered 
situation, now in presence of a GSL with finite strength 
V = TT, with the same lattice geometry and energy as used 
for Fig. 3, but relatively strong amplitude fluctuations of 
±30%. A strong filtering effect analogous to panel b is 
observed, where the largest part of the transmitted prob- 
ability density is concentrated in the forward direction 
9 = 0, as expressed by the overall tan^ 0-behavior of the 
Lyapunov exponent, Eq. (40). 

The wave propagation in the clean GSL of panel c is 
quite isotropic, because the associated dispersion relation 
is almost isotropic for the parameters chosen. If, however, 
the potential strength of a scalar GSL is greater than a 
certain critical value, new Dirac points emerge.""" The 
resulting, strongly anisotropic dispersion relation then 
collimates the wave packet."' This is shown in panel e, 
where the potential strength v = At: makes the wave 
packet stay sharply focused in the forward direction. We 
have investigated whether this collimation effect is robust 
against disorder. Panel f shows the effect of 10% fluc- 
tuations in potential strength. The collimation is seen 
to survive, with hardly noticeable disorder effects. A 
quantitative analysis is difficult because the new Dirac 
points appear at band edges where the weak-disorder ex- 
pansion we have used fails. Instead, one could possibly 
adapt the appropriate singular-point expansions''"''"'" to 
the Dirac-GSL problem, which is a research program be- 
yond the scope of the present work. Here, we conclude 
that disorder-induced filtering can coexist with band- 
structure collimation. 
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VI. CONCLUDING REMARKS 

Drawing on a general weak-disorder expansion, we 
have derived the Lyapunov exponent (inverse localiza- 
tion length) of various ID disordered GSLs modeled by 
random delta or rectangular potentials, both for scalar 
and vector potentials. The analytical results have been 
thoroughly checked by numerical experiments. We em- 
phasize that, though the GSL is assumed to be ID, the 
physics is far more complicated than for a conventional 
ID scattering problem due to the intrinsic coupling be- 
tween the translational motion and the spinor degree of 
freedom. One important complication we have predicted 
is the strong dependence of the localization length on the 
incident angle of the charge carriers injected to a GSL. 
To our knowledge, this is the first time that a complete 
theoretical picture of this incident-angle dependence is 
obtained. We have also proposed to exploit such angu- 
lar dependence of the localization length to turn disorder 
into good use, namely, a possible disorder-assisted fil- 
tering effect. Considering that large-size GSLs may be 
manufactured in the near future, our theoretical results 
offer a quantitative tool to analyze and predict disorder 
effects in GSLs. 

Our analytical and numerical results also provide ev- 
idence for intriguing delocalization resonances: Along 
specific incident angles, the localization exponent can be 
identically zero, or at least approach zero for weak dis- 
order. Both scalar and vector GSLs admit delocalization 
resonances in the conduction band, but for opposite rea- 
sons: scalar potentials can have an approximate, weak- 
disorder resonance because a complex term has zero so- 
lutions [i.e., 5" = 0, see Eq. (40) and (41)], whereas vec- 
tor potentials have an approximate resonance because 
of a simple barrier resonance condition [sint/j = 0, see 
Eq. (55)]. Moreover, the corresponding (5-limits of scalar 
and vector GSLs show very distinct features: the scalar 
(5GSL admits an exact delocalization by virtue of an inter- 
peak resonance, whereas the vector (5GSL has no reso- 
nance at all in the conduction band. In all cases, it is 
important to realize that whenever numerical or labo- 
ratory experiments are performed with finite-size sam- 
ples, a lowest-order vanishing Lyapunov exponent can 
very well appear as a rather sharp mobility jump, which 
signals an effective delocalization across the sample.'""'''' 

In the context of 2D GSLs, a recent study'" cautioned 
that lattice constants less than 10 nm may induce inter- 
valley scattering or sublattice symmetry breaking, either 
of which may lead to a band gap and hence break the 
linear dispersion relation of the charge carriers. The im- 
plication of this important finding for our work is twofold. 
First, to directly apply our theoretical results based on 
a linear dispersion relation, it is safer to consider GSLs 
with lattice constants larger than 10 nm or with a po- 
tential preserving the symmetry between different Dirac 
points or between different sublattices. Second, as a pos- 
sible extension of this work, one may now also apply our 
main theoretical tool here to investigate how a disordered 



GSL with a sufficiently small lattice constant may gen- 
erate a novel physical situation, where charge carriers 
possess disordered mass as a consequence of inter-valley 
scattering or sublattice symmetry breaking. 
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Appendix A: Details of weak-disorder expansion 

This appendix provides some details of the analytical 
calculation leading to the weak-disorder Lyapunov expo- 
nent given by Eq. (22) and Eq. (23). 

1. Absence of mixed- fluctuation terms 



The starting point is Eq. (21), where X± 
R, describes a propagating solution in the clean GSL. A 
Taylor expansion to quadratic order in the fluctuations 
e„ leads to 



1 

2N 



ln|(Pw)i 



=Re. 



N 



Re- 



(Al) 



2(n— m) 



First, we justify that the last line only gives a negligible 
contribution under the ensemble average. In terms of 
the complex random variable z„ = fnA^" — e„e^™'^, the 
double sum rewrites 



n<m 



1 

m 



N 



1 

m 



N 



In the second term, we recognize the variance \z\ 

in the limit TV — cx). The whole expression (A2) can 
be written as (T^(|?/ip — l)/2, where the random variable 

Ui = J2n=i -^ra V(V^ cr) fluctuates as samples i are drawn 
from the ensemble. Now, according to the Berry-Esseen 
theorem, in the limit N ^ oo the probability distribution 
of |y| converges to the standard normal distribution, with 
unit variance — 1. As a consequence, cr^dyp — 1) = 0, 
such that the whole expression (A2) gives zero contribu- 
tion after the ensemble average. 

Then, the vanishing fluctuation mean (19) makes also 
the first line in (Al) vanish. Thus, only the variance (20) 
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in the second line contributes to 



Elementary algebra leads to 



7 — lim 



ln|(Pjv)ii|^ 
2N 



-Re 



which is the result stated as Eq. (22). 



(A3) 



2. Diagonalization procedure 

As the last task, we need to express the matrix el- 
ements M'li and of the diagonal representation 
through the transfer-matrix parameters {a,/3, 0} as de- 
fined in (12). In an intermediate step, we parameterize 
the transfer matrix as 



M ^ 



where 



b = e'P tan( 



(A4) 



(A5) 
(A6) 



have to satisfy the constraint det Af — |ap — |&p = 1. The 
diagonal representation M = diag(A+, A_) = P^^MP is 
attained by a basis transformation with 



P 



h b 

— a A_ — a 



1 



A_ — a 



fe(A_ - A+) \-{\+-a) b 



(A7) 
(A8) 



Again, we assume that the eigenvalues A± = e^'^ form a 
complex conjugate (/j, € R), and non-degenerate (/x ^ nir 
for all n S S) pair since we seek the Lyapunov exponent 
of inside-conduction-band solutions. 



= {p-Hrpu = ?^±M^, (A9) 

where the useful identities A+ + A_ = a + a* and 
(detM)' = 2Re{a*a' - b*b'} = have been employed. 
Furthermore, Re{a} ~ cos/i entails Re{a'} — — /i'sin/i, 
and since A+ — A_ = 2zsin/i, Eq. (A9) implies 



(a/;,/Ah 



sin /i 



which is the second term needed in Eq. (A3). 
Proceeding similarly, one finds for the first term 



M('i _ Re{a"} Re{c 



h*h"} 



A+ i sin /i — isin/i 

We only need its real part, 



(cos /i — i sin /i). 

(All) 



Re |a;^M{'i} = Re{a*a" - b*b") = \b' 

= (0'2 - a'2) sec^ (j) + P'^ tan^ 0. (A12) 
Substituting Eq. (AIO) and (A12) into Eq. (A3), we have 



^ f (2 , / 1,12 I2\ 



a'^) sec^ (/) + /?'^ tan^ <p} , (A13) 
Further algebraic manipulations lead to the identity 



,'2 



J,'2 „'2 



a'^)sec^ 



tan 



sm a 
am 6 



1 — sec^ cos^ a 

(A14) 

which then results in the final expression Eq. (23) for the 
weak-disorder Lyapunov exponent. 
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